Workshop4_Spatial

# Loading packages

#library(tidyverse)
#install.packages("sf") 
#install.packages("terra")
#install.packages("tmap")
library(tidyr)
library(sf) # simple features

## Linking to GEOS 3.11.0, GDAL 3.5.3, PROJ 9.1.0; sf_use_s2() is TRUE

library (terra) # for raster

## terra 1.7.71

## 
## Attaching package: 'terra'

## The following object is masked from 'package:tidyr':
## 
##     extract

library(tmap) # Thematic maps are geographical maps in which spatial data distributions are visualized

## The legacy packages maptools, rgdal, and rgeos, underpinning the sp package,
## which was just loaded, will retire in October 2023.
## Please refer to R-spatial evolution reports for details, especially
## https://r-spatial.org/r/2023/05/15/evolution4.html.
## It may be desirable to make the sf package available;
## package maintainers should consider adding sf to Suggests:.
## The sp package is now running under evolution status 2
##      (status 2 uses the sf package in place of rgdal)

## Breaking News: tmap 3.x is retiring. Please test v4, e.g. with
## remotes::install_github('r-tmap/tmap')

library(readr)
dat <- read_csv("../data/copepods_raw.csv")

## Rows: 5313 Columns: 11
## ── Column specification ────────────────────────────────────────────────────────
## Delimiter: ","
## chr (5): sample_time_utc, project, route, vessel, region
## dbl (6): silk_id, segment_no, latitude, longitude, meanlong, richness_raw
## 
## ℹ Use `spec()` to retrieve the full column specification for this data.
## ℹ Specify the column types or set `show_col_types = FALSE` to quiet this message.

# Creating a figure that plots our spatial data, but isnt neccessarily spatially accurate
library(ggplot2)
ggplot(dat) + 
  aes(x = longitude, y = latitude, color = richness_raw) +
  geom_point()

#So, now let’s look at the richness data 

ggplot(dat, aes(x = latitude, y = richness_raw)) + 
  stat_smooth() + 
  geom_point()

## `geom_smooth()` using method = 'gam' and formula = 'y ~ s(x, bs = "cs")'

#turn our point data into a spatially referenced data frame

sdat <- st_as_sf(dat, coords = c("longitude", "latitude"), 
                 crs = 4326)

#Coordinate reference systems
#Our copepod coordinates are long-lat, so we chose a common ‘one-size-fits-all’ GCS called WGS84 to define the crs using the EPSG code 4326.

crs4326 <- st_crs(4326)
crs4326 # look at the whole CRS

## Coordinate Reference System:
##   User input: EPSG:4326 
##   wkt:
## GEOGCRS["WGS 84",
##     ENSEMBLE["World Geodetic System 1984 ensemble",
##         MEMBER["World Geodetic System 1984 (Transit)"],
##         MEMBER["World Geodetic System 1984 (G730)"],
##         MEMBER["World Geodetic System 1984 (G873)"],
##         MEMBER["World Geodetic System 1984 (G1150)"],
##         MEMBER["World Geodetic System 1984 (G1674)"],
##         MEMBER["World Geodetic System 1984 (G1762)"],
##         MEMBER["World Geodetic System 1984 (G2139)"],
##         ELLIPSOID["WGS 84",6378137,298.257223563,
##             LENGTHUNIT["metre",1]],
##         ENSEMBLEACCURACY[2.0]],
##     PRIMEM["Greenwich",0,
##         ANGLEUNIT["degree",0.0174532925199433]],
##     CS[ellipsoidal,2],
##         AXIS["geodetic latitude (Lat)",north,
##             ORDER[1],
##             ANGLEUNIT["degree",0.0174532925199433]],
##         AXIS["geodetic longitude (Lon)",east,
##             ORDER[2],
##             ANGLEUNIT["degree",0.0174532925199433]],
##     USAGE[
##         SCOPE["Horizontal component of 3D system."],
##         AREA["World."],
##         BBOX[-90,-180,90,180]],
##     ID["EPSG",4326]]

crs4326$Name # pull out just the name of the crs

## [1] "WGS 84"

crs4326$wkt  # crs in well-known text format

## [1] "GEOGCRS[\"WGS 84\",\n    ENSEMBLE[\"World Geodetic System 1984 ensemble\",\n        MEMBER[\"World Geodetic System 1984 (Transit)\"],\n        MEMBER[\"World Geodetic System 1984 (G730)\"],\n        MEMBER[\"World Geodetic System 1984 (G873)\"],\n        MEMBER[\"World Geodetic System 1984 (G1150)\"],\n        MEMBER[\"World Geodetic System 1984 (G1674)\"],\n        MEMBER[\"World Geodetic System 1984 (G1762)\"],\n        MEMBER[\"World Geodetic System 1984 (G2139)\"],\n        ELLIPSOID[\"WGS 84\",6378137,298.257223563,\n            LENGTHUNIT[\"metre\",1]],\n        ENSEMBLEACCURACY[2.0]],\n    PRIMEM[\"Greenwich\",0,\n        ANGLEUNIT[\"degree\",0.0174532925199433]],\n    CS[ellipsoidal,2],\n        AXIS[\"geodetic latitude (Lat)\",north,\n            ORDER[1],\n            ANGLEUNIT[\"degree\",0.0174532925199433]],\n        AXIS[\"geodetic longitude (Lon)\",east,\n            ORDER[2],\n            ANGLEUNIT[\"degree\",0.0174532925199433]],\n    USAGE[\n        SCOPE[\"Horizontal component of 3D system.\"],\n        AREA[\"World.\"],\n        BBOX[-90,-180,90,180]],\n    ID[\"EPSG\",4326]]"

#Feature collection (points)
sdat

## Simple feature collection with 5313 features and 9 fields
## Geometry type: POINT
## Dimension:     XY
## Bounding box:  xmin: 89.6107 ymin: -65.2428 xmax: 174.335 ymax: -16.80253
## Geodetic CRS:  WGS 84
## # A tibble: 5,313 × 10
##    silk_id segment_no sample_time_utc  project route vessel      meanlong region
##  *   <dbl>      <dbl> <chr>            <chr>   <chr> <chr>          <dbl> <chr> 
##  1       1          1 26/06/2009 22:08 AusCPR  BRSY  ANL Windar…     153. East  
##  2       1          5 26/06/2009 23:12 AusCPR  BRSY  ANL Windar…     153. East  
##  3       1          9 27/06/2009 0:17  AusCPR  BRSY  ANL Windar…     153. East  
##  4       1         13 27/06/2009 1:22  AusCPR  BRSY  ANL Windar…     153. East  
##  5       1         17 27/06/2009 2:26  AusCPR  BRSY  ANL Windar…     153. East  
##  6       1         18 27/06/2009 2:43  AusCPR  BRSY  ANL Windar…     153. East  
##  7       1         26 27/06/2009 4:52  AusCPR  BRSY  ANL Windar…     153. East  
##  8       1         30 27/06/2009 5:57  AusCPR  BRSY  ANL Windar…     153. East  
##  9       1         33 27/06/2009 6:45  AusCPR  BRSY  ANL Windar…     153. East  
## 10       1         37 27/06/2009 7:50  AusCPR  BRSY  ANL Windar…     153. East  
## # ℹ 5,303 more rows
## # ℹ 2 more variables: richness_raw <dbl>, geometry <POINT [°]>

# Cartography

# Plotting Richness
plot(sdat["richness_raw"])

plot(sdat)

#Thematic maps for communication

#Build and add on layers using using tmap 

tm1<- tm_shape(sdat) + 
  tm_dots(col = "richness_raw")

# Saving our figure
tmap_save(tm1, filename = "Richness-map.png", 
          width = 600, height = 600)

## Map saved to /Users/simonelinderman/Desktop/Marine Tech 1/Data Science/Module4_Personal/code/Richness-map.png

## Resolution: 600 by 600 pixels

## Size: 2 by 2 inches (300 dpi)

# Mapping spatial polygons as layers

# First, load spatial data
aus <- st_read("../data/spatial-data/Aussie/Aussie.shp")

## Reading layer `Aussie' from data source 
##   `/Users/simonelinderman/Desktop/Marine Tech 1/Data Science/Module4_Personal/data/spatial-data/Aussie/Aussie.shp' 
##   using driver `ESRI Shapefile'
## Simple feature collection with 8 features and 1 field
## Geometry type: MULTIPOLYGON
## Dimension:     XY
## Bounding box:  xmin: 112.9211 ymin: -43.63192 xmax: 153.6389 ymax: -9.229614
## Geodetic CRS:  WGS 84

shelf <- st_read("../data/spatial-data/aus_shelf/aus_shelf.shp")

## Reading layer `aus_shelf' from data source 
##   `/Users/simonelinderman/Desktop/Marine Tech 1/Data Science/Module4_Personal/data/spatial-data/aus_shelf/aus_shelf.shp' 
##   using driver `ESRI Shapefile'
## Simple feature collection with 1 feature and 1 field
## Geometry type: MULTIPOLYGON
## Dimension:     XY
## Bounding box:  xmin: 112.2242 ymin: -44.1284 xmax: 153.8942 ymax: -8.8798
## Geodetic CRS:  GRS 1980(IUGG, 1980)

# Mapping your polygons

tm_shape(shelf) + 
  tm_polygons()

# Layering copepod data points
tm_shape(shelf, bbox = sdat) + 
  tm_polygons() +
  tm_shape(aus) + 
  tm_polygons() + 
  tm_shape(sdat) + 
  tm_dots()

#Exploring t_map
vignette('tmap-getstarted')

## starting httpd help server ... done

tm_shape(shelf) + 
  tm_polygons()

# Adding additional element, such as colors
tm_shape(shelf, bbox = sdat) + 
  tm_polygons() +
  tm_shape(aus) + 
  tm_polygons() + 
  tm_shape(sdat) + 
  tm_dots(size = 0.01, col = "blue")

# Vignette Exercises

# World map 
data("World") +
  tm_shape(World) +
    tm_polygons("HPI")

# Multiple shapes and layers

data(World, metro, rivers, land)

tmap_mode("plot")

## tmap mode set to plotting

## tmap mode set to plotting
tm_shape(land) +
    tm_raster("elevation", palette = terrain.colors(10)) +
tm_shape(World) +
    tm_borders("white", lwd = .5) +
    tm_text("iso_a3", size = "AREA") +
tm_shape(metro) +
    tm_symbols(col = "red", size = "pop2020", scale = .5) +
tm_legend(show = FALSE)

## old-style crs object detected; please recreate object with a recent sf::st_crs()

## old-style crs object detected; please recreate object with a recent sf::st_crs()
## old-style crs object detected; please recreate object with a recent sf::st_crs()

# Facets
tmap_mode("view")

## tmap mode set to interactive viewing

tm_shape(World) +
    tm_polygons(c("HPI", "economy")) +
    tm_facets(sync = TRUE, ncol = 2)

#Splitting the spatial data
tmap_mode("plot")

## tmap mode set to plotting

## tmap mode set to plotting

data(NLD_muni)
tm1 <- tm_shape(NLD_muni) + tm_polygons("population", convert2density = TRUE)
tm2 <- tm_shape(NLD_muni) + tm_bubbles(size = "population")

tmap_arrange(tm1, tm2)

## old-style crs object detected; please recreate object with a recent sf::st_crs()

## old-style crs object detected; please recreate object with a recent sf::st_crs()
## old-style crs object detected; please recreate object with a recent sf::st_crs()
## old-style crs object detected; please recreate object with a recent sf::st_crs()
## old-style crs object detected; please recreate object with a recent sf::st_crs()
## old-style crs object detected; please recreate object with a recent sf::st_crs()
## old-style crs object detected; please recreate object with a recent sf::st_crs()
## old-style crs object detected; please recreate object with a recent sf::st_crs()
## old-style crs object detected; please recreate object with a recent sf::st_crs()
## old-style crs object detected; please recreate object with a recent sf::st_crs()
## old-style crs object detected; please recreate object with a recent sf::st_crs()
## old-style crs object detected; please recreate object with a recent sf::st_crs()
## old-style crs object detected; please recreate object with a recent sf::st_crs()
## old-style crs object detected; please recreate object with a recent sf::st_crs()
## old-style crs object detected; please recreate object with a recent sf::st_crs()
## old-style crs object detected; please recreate object with a recent sf::st_crs()
## old-style crs object detected; please recreate object with a recent sf::st_crs()
## old-style crs object detected; please recreate object with a recent sf::st_crs()
## old-style crs object detected; please recreate object with a recent sf::st_crs()
## old-style crs object detected; please recreate object with a recent sf::st_crs()

## Legend labels were too wide. Therefore, legend.text.size has been set to 0.32. Increase legend.width (argument of tm_layout) to make the legend wider and therefore the labels larger.

#Basemaps and overlay tile maps

tmap_mode("view")

## tmap mode set to interactive viewing

tm_basemap("Stamen.Watercolor") +
tm_shape(metro) + tm_bubbles(size = "pop2020", col = "red") +
tm_tiles("Stamen.TonerLabels")

## old-style crs object detected; please recreate object with a recent sf::st_crs()

## old-style crs object detected; please recreate object with a recent sf::st_crs()
## old-style crs object detected; please recreate object with a recent sf::st_crs()
## old-style crs object detected; please recreate object with a recent sf::st_crs()
## old-style crs object detected; please recreate object with a recent sf::st_crs()

## Legend for symbol sizes not available in view mode.